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We discuss computational analysis of phononic thermal conduction in the suspended 
membrane geometry, in the experimentally commonly appearing case where heat can 
flow out radially in two dimensions from a central source. As we are mostly interested 
in the low-temperature behavior where bulk scattering of phonons becomes irrelevant, 
we study the limit where all phonon scattering takes place at the membrane surfaces. 
Moreover, we limit the discussion here to the case where this surface scattering is fully 
diffusive, the so called Casimir limit. Our analysis shows that in the two-dimensional 
case, no analytic results are available, in contrast to the well known ID Casimir limit. 
Numerical solutions are presented for the temperature profiles in the membrane radial 
direction, for several different membrane thicknesses. Our results can be applied, for 
example, in the design of membrane-supported bolometric radiation detectors. 
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I. INTRODUCTION 



Thermal conduction in insulators takes place via phonon transport, which is typically 
limited in bulk by various scattering mechanisms, with phonon-phonon scattering being 
dominant in the room temperature region^. In this standard case, the thermal conductance 
G, defined as the coefficient between conducted power P and temperature difference, P = 
GAT, can be written in terms of the local thermal conductivity k and geometry dependent 
terms. However, if one goes to temperatures much below the Debye temperature T < 9 D /30, 
the phonon-phonon scattering dies away 1 , and in pure enough samples one is left with just 
the scattering at the boundary of the sample. In that case, therefore, there is no true local 
thermal conductivity left. Nevertheless, data is still often plotted in those units, and one 
has to remember that the thermal conductivity actually depends on the sample dimensions 
in that case. 

The boundary scattering limit was first analyzed for the case where heat is flowing along 
a one-dimensional rod log time ago by Casimir-. In his analysis, the assumption was made 
that the scattering was fully random at the surface, in other words a phonon impinging 
on the surface scatters diffusively to any direction in the half-space with equal probability. 
He obtained the simple and intuitive result that for an infinitely long cylindirical rod, the 
thermal conductance can be written in terms of the same formula as for bulk scattering, 
k = ^Cvl, where C is the specific heat capacity of the phonon gas, v is a properly defined 
average speed of sound 3 , but with the mean free path I simply given by the diameter of 
the rod. This diffusive boundary scattering case was later dubbed the Casimir limit, to 
differentiate it from more general cases, where the boundary scattering is partially or fully 
specular-^. 

In many experimental cases the heat flow is of course one-dimensional, and the standard 
Casimir results can be compared for example to results on phononic thermal conduction 
in one-dimensional nanowires at low temperatures^^ - —. However, in some cases the ex- 
perimental realization is two-dimensional for the heat flow. By two-dimensional we here 
mean that there is a central heat source, and heat can spread uniformly in all radial di- 
rections. This is particularly relevant for the case of membrane-supported low-temperature 
bolometric radiation detectors^, where the phononic thermal conductance of a suspended 
membrane is the limiting mechanism for heat dissipation. Recent experiments on suspended 
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SiN membranes^ have clearly demonstrated that the sub-Kelvin thermal conductance in the 
membrane geometry is not limited by bulk scattering, by studying the temperature profile 
in the membrane directly This also agrees with earlier results where different samples were 
compared 1 -, or where the absolute value of thermal conductance was seen to approach the 
ballistic limit (fully specular surface scattering)^. For full understanding of the 2D experi- 
ments, one certainly should also understand the Casimir limit for 2D heat flow. Somewhat 
surprisingly, we did not find any discussion of it in literature. 

In this paper, we thus answer the question how the 2D heat flow behaves in the Casimir 
limit from a theoretical point of view. We present only a numerical solution: the equation 
that we derive does not allow an analytical solution, even in the simplest, most symmetric 
geometry. It remains to be seen whether an analytical solution exists, but based on the 
discussion here it seems unlikely. There is therefore a major difference between the 2D and 
ID Casimir results, and most importantly, one cannot use the bulk formulas for thermal 
conductivity at all in the 2D case. In other words, there is no equivalent general rule that 
the mean free path would be proportional to some dimension of the sample. 

II. THERMAL CONDUCTANCE LIMITED BY BULK SCATTERING IN 
2D 

Before we present the equations for 2D Casimir-limited phonon transport, let us discuss 
the opposite case of 2D heat flow limited by bulk scattering^^. In ID everyone is familiar 
with how the extensive variable thermal conductance G depends on the intensive thermal 
conductivity k: G = kA/L, where A is the cross sectional area and L is the length of the 
sample. What about in 2D? In that case heat flows out radially, and the cross-sectional 
area is a function of the radial coordinate. Clearly, the equation between G and k is more 
complicated. We can derive it by using the Fourier law for the heat flux j = — VT, and since 
we assume cylindrical symmetry, only the radial derivative in the gradient of the temperature 
remains, so that the power is given by 



By noting that A(r) = 2nrd, where d is the membrane thickness, separating the variables 
r and T, integrating, and assuming small temperature difference AT << T, one can write 




(1) 
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x cLS cL function of AT, and thus identify G: 



G 



2nd 



(2) 



ln(ri/r ) 



where r± > r are the two radii where the two temperatures in AT are defined. One sees 
that there is an explicit dependence on the two coordinates left: G is not a function of just 
the difference of r% and T2 like in the ID case. If, furthermore, we assume that we are far 
from the origin t\ — r << r , we can Taylor expand the logarithm and simplify Eq. [2] to 
get 



where L = r 1 — r . This is finally equivalent with the ID A = 27rr d. We 

thus see what we intuitively expect: for small L far from origin, the curvature difference 
between the two circles is small so that the 2D case approaches the ID case. However, in 
many experimental situations this limit is not reached and the full formula (Eq. [2]) must 
be used to determine k from the measured G. Moreover, when the temperature difference 
is not small, one cannot define a linear thermal conductance in terms of the temperature 
difference anymore. In that case, it is more useful to calculate the power P, which is an 
explicit function of both the local temperature and the bath temperature T&. One can still 
define a differential thermal conductance Gd = dP/dT which agrees with the result in Eq. 
|2j If one assumes some functional form for the temperature dependence of the thermal 
conductivity, for example k = aT m , one can invert the power vs. T relation, arriving^- 2 at a 
simple result for the temperature profile in the bulk scattering case: 



where T(ri) is a know temperature at radius r±, for example at the membrane edge where 
it would equal the bath temperature. 

III. THE APPROACH TO STUDY THE CASIMIR LIMIT IN 2D 

As in Casimir's original work 2 , we also take the same approach to model a fully diffusive 
surface. When the scattering on the surfaces is fully diffusive, the phonon emission from the 
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surface element follows the same laws as photon blackbody radiation, but with a redefined 
phononic Stefan-Bolzmann constant a. For an isotropic material, a is given 1 ^ by 



where c t and q are the transverse and longitudinal speeds of sound of the material, 
respectively. These are the only material parameters in the problem. From the analogy 
with electromagnetic radiation 1 ^, we therefore know that the total hemispherically emitted 
power of a surface element is cxT 4 , and that the emitted intensity from a surface element dS 
into solid angle dfl is given by Lambert's law: dl = ^T 4 cosOdSdQ, where 9 is the angle that 
the directional vector makes with the surface normal. In steady state, there must exist a 
radiative balance between the emitted and absorbed power for each surface element, leading 
to an equation ^TfdSi J2j cos9idQj = ^ J2j Tf cos OjdSjdQi + qdSi, where the left side is the 
total emitted power of element dSi into surface elements dSj, the first term on the right side 
is the total absorbed power at element dSi from all other surface elements dSj, and the last 
term is an external power absorption due to a metallic heater, for example (q has units of 
W/m 2 ). The solid angle terms can be written in terms of the two surface elements and are 
thus dVLj = dSj cos Oj/Rfj and dfli = dSi cos 9i/Rfj, where Rij is the distance between the 
two surface elements. 

As we are interested in the 2D case with cylindrical symmetry, we can write the above 
equation in cylindrical coordinates and go to the continuum limit by converting the sums to 
integrals. At this point, we should also explicitly take into account of the membrane geom- 
etry. When considering a surface element on the membrane, we can see that there are two 
different contributions to the dSj elements: one is the opposing, parallel membrane surface 
(no power is exchanged between the surface elements on the same side due to Lambert's 
law), and the other is the membrane edge at radial coordinate R, which has a different con- 
tribution because the edge surface is perpendicular to the membrane surfaces. To maintain 
the cylindrical symmetry, we have to assume a circular membrane edge and a circular heater 
geometry. Moreover, we simplify the discussion slightly by also assuming that the heater is 
symmetric with respect to the direction perpendicular to the membrane (i.e both top and 
bottom surfaces are heated): in this case the upper and lower membrane surfaces have iden- 
tical temperatures everywhere. This assumption could actually be relaxed by writing the 
radiative balance equations for the top and bottom surfaces separately, which would lead to 




(5) 
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a coupled system of equations, which at this point is an unnecessary complication. We also 
take the temperature of the membrane edge to be fixed by the heat bath, as it is in direct 
contact with the bulk of the sample. This is a good approximation when the membrane 
thickness d is much larger than the dominant wavelengths of the thermal phonons involved; 
for very thin membranes there should be a corrections due to diffractive processes. 

Using the above assumptions, we arrive at an integral equation desribing the radiative 
balance for 2D Casimir limit as a function of the radial coordinate r\ : 

^ rR sin 2 # a (** r d sin 6 b (R 2 - r x R cos (p 2 )\ 



m 



b 



T 4 (ri) [ f d<p 2 f dr 2 r 2 Sm a + [ dfc f dz 
\Jo Jo R a Jo Jo 

= d<p 2 dr 2 T (r 2 )r 2 p2 + T (R) / d<f> 2 / dz -3 + ~qf(r x y$) 

Jo Jo Jo Jo Rl a 

where the first term on the left is the emitted power from one parallel membrane surface 
to another and the second from the membrane surface to the membrane edge at r\ = R. 
On the right we have the absorption terms, the first from the other parallel surface, the 
second from the membrane edge and the third the direct power input, which now has a 
radial profile f{r\). Here 9i are the angles between the membrane surface and the distance 
vectors (with lengths R a and Rf,) in the two cases. From geometry we can easily derive that 
R 2 =r 2 + r 2_ 2ri r 2 cos(0i - 2 ) + d 2 , sin 6 a = d/R a , R 2 = r 2 + R 2 - 2r x R cos(0i - 2 ) + z 2 , 
and sin^t = z/R^. Thus, it becomes a matter of elementary integration to simplify Eq. [6l 
resulting in the following final linear integral equation for the temperature profile Z(r\) = 
T\n): 

Z( n ) = f R dr 2 G( ri ,r 2 )Z(r 2 ) + Z(R)H(n) + Cf{r x ), (7) 
Jo 

where the kernel G(ri,r 2 ) is given by 

G(r 1 ,r 2 ) = ^ 2 r 2 (r 2 + r 2 + d 2 ) 

\{r 2 + r 2 + d 2 ) 2 -Ar 2 r 2 f 12 

the "edge" function H(r\) as 

/ n 1 I r 2 + d 2 -R 2 \ 

H(ri) = - 1 + 1 , (9) 

2 \^r 2 + R 2 + d 2 ) 2 - ArfR 2 J 

and the constant C is the normalized external power input C = 2q/a so that the function 
f(r±) is one where power is applied and zero where not. 
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IV. NUMERICAL RESULTS 



We choose to solve equation [7| numerically, as no obvious analytic solution exists. Using 
the terminology of mathematical literature, equation [7] can be classified as a linear Fredholm 
equation of the second kind^. For the numerics it is very important how the kernel behaves. 
As is seen from Eq. |H1 it is not symmetric, nor is it a function of just a difference of the 
variables. It is, however, non-singular at n = r2, so that standard numerical techniques 
for these type of integral equations work. We have chosen to use the so called Nystrom 
method^, which uses the Gauss- Legendre quadrature rule for discretization of the integral, 
and trianglular decomposition techniques for the inversion of the obtained linear equations. 
By trial and error we have seen that the details of the problem influence how many points 
are needed in the discretization. The higher the aspect ratio between membrane radius and 
thickness, the higher number of points must be used. Here, we we had to go up to 900 points 
(900 by 900 matrix inversion) to get accurate results for the thinnest membranes studied 
here (aspect ratio R/d = 250). We note that thinner membranes would require even more 
points. 

In Fig. [1] we present the results for the temperature profiles in the radial direction for 
several different aspect ratios ranging from R/d = 2.5 to R/d = 250, keeping R constant, in 
linear (a) and log-log scales (b), with a heater radius = 1. As there is no intrinsic length 
scale in the problem, all results can be plotted in arbitrary units; in other words the plots 
could represent results for R = 250/mi and d ranging from 1 /zm to 100 /xm, or any other 
length scale. Temperature is also plotted in the naturally scaled units, which depend on the 
material through a and heater power through q, and heater powers are selected in such a 
way that the temperature at the center r = is the same for all membrane thicknesses. In 
addition to the numerical Casimir-limit results, we also plot the bulk diffusive case of Eq. 
® 

First, we see that the bulk diffusive result gives higher temperatures than all the calculated 
Casimir cases (except very near the mebrane edge at r = 250), and that the Casimir results 
clearly produce different shapes for the profiles. The profiles also clearly depend on the 
aspect ratio, so that there is no size-independent universal shape for the Casimir case. This 
is in contrast with the bulk diffusive and ballistic limits 12 . The second very obvious feature 
is that the temperature profile is not necessarily continuous, and there are jumps both at 
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FIG. 1. [Color online] Calculated temperature profiles for 2D Casimir heat conduction (points) 
for R = 250 and d = 1,3.16, 10,31.6, 100, where R and d have arbitrary units in (a) linear scale, 
(b) log-log scale. Higher curves have lower d. Dashed line shows the bulk diffusive result, Eq. [H 
and dotted line the value of bath temperature used. In (b), the solid lines are power law fits with 
exponents -1/9, -1/8, -1/6, -1/5 from top to bottom. 

the heater edge at r = 1 and at the membrane edge at r = 250. This is in agreement with 
Monte Carlo simulations of the 1-D Casimir limits, where temperature jumps at the rod 
ends were also seen. The temperature jump at the heater edge grows significantly in size 
when the membrane becomes thicker, so much so that in the thickest membranes studied 
here, a whopping 93 % of the temperature drop in the sample occurs right at the heater 
edge. The drop at the membrane edge is much more modest in all the cases, now being 
largest ( 9 %) for the thinnest membrane. The growth of the drop at the heater edge is 
intuitively understandable: for thick membranes the heater can emit to a much larger solid 
angle, therefore distributing the heat load much more widely and non-locally, resulting in 
much cooler membrane outside of heater area. 

By plotting the profiles in log-log scale, Fig. Q](b), one can also try to infer phenomenolog- 
ical dependecies for the temperature decay as a function of radial distance. The Casimir-limit 
temperature profiles outside the heater area seem to have several regions, depending on d: 
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FIG. 2. [Color online] Calculated temperature profiles for 2D Casimir heat conduction (points) for 
R = 250, d = 1, and r^=l, with varying heat input power density 2q/a = 10 -4 , 10~ 3 , 0.01, and 
0.1 K 4 in (a) linear scale, (b) log- log scale. Higher curves have higher q. The dotted line is the 
value of bath temperature 0.1 K used. 

there is an initial flat region for the thicker membranes (extending up to r ~ 30 for d = 100), 
then an intermediate power law decay, where the exponent is very weak and depends on d, 
and then finally a stronger decay beyond r ~ 50, resembling the bulk diffusive limit. We 
have fitted simple power laws r~ b to some of the intermediate regions of the plots, where 
b was varying from 1/9 to 1/5. However, no clear power-law region exisits for the thickest 
membrane result. 

The profiles presented in Fig. [1] are universal in terms of the power input (as implied by 
the scaling used for T) as long as the temperatures at the membrane edge remain higher 
than the bath temperature T^ dge » T^ ath . If that is not the case, the shapes are distorted 
and become flatter. We show this effect in Fig. |2j where results for d = 1 with low input 
power levels are shown. Clearly, with low input power, the temperature profile is very flat 
and does not show a discontinuous drop at the membrane edge, until a higher power level 
is reached (2g/er = 0.1 for this aspect ratio). The power law decay, on the other hand, is 
not strongly affected until at the very lowest power input levels. The strongest difference 
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FIG. 3. [Color online] Calculated temperature profiles for 2D Casimir heat conduction (points) 
for R = 250, d = 1, and constant total heating power r^lqja = 490 K 4 (arb) 2 , but with varying 
heater radius =1,3,7,25,50,100,150 in (a) linear scale, (b) log-log scale. The dashed lines mark 
the values of the heater radii. 

between the low power profiles (Fig. [2]) and the high power profile (Fig. [1]) is in the third 
region beyond r = 50, where the stronger, bulk-like temperature decay disappears for low 
input powers, and becomes linear. Similar behavior takes place for thicker membranes (not 
shown) . 

Finally, we want to investigate how the temperature profiles depend on the heater size. 
Fig. [3] shows the calculated temperature profiles for d = 1, R = 250, with varying heater 
size from — 1 to — 150, but by keeping the total power constant. Clearly, all the curves 
come together outside of the heater region; only the temperature profile within the heater are 
is affected, with heater temperature naturally rising with rising power density. We also point 
out that the temperature wihin the heater region is in no way constrained to be constant. In 
fact, as the larger heater size results clearly show, there is a sizeable temperature gradient 
wihin the heater area. The geometry only demands that the temperature gradient at the 
center r = is zero. Also, the finite temperature drop at the heater edge disappears for 
larger heater radii r h > 3. 
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V. DISCUSSION AND CONCLUSIONS 



We studied theoretically two-dimensional phonon transport in suspended membranes 
in the limit where all the scattering takes place on the membrane surfaces, and with the 
assumption that all the scattering is diffuse, the so called Casimir limit. What we mean by 
two-dimensionality in this context is that there is a central heat source, from which heat 
emanates, spreading out into the membrane and the heat bath, which is located at the edge 
of the membrane some distance away. We showed that already in the more general case of 
2D heat low with bulk scattering, there is no way to describe in general the heat transport 
properties with just thermal conductivity, as thermal conductance depends explicity on the 
positions where temperature is measured. We also showed that in the Casimir limit, the 
situation is complicated even further so that there is no universal temperature profile, but 
instead it depends on the aspect ratio between the membrane radius and thickness. We 
presented the derivation of an integral equation for the temperature profile from power 
balance at the surfaces, and proceeded to solve the equation numerically for several cases. 

Our formulation assumed a full cylindrical symmetry for simplification. However, in 
reality suspended membranes tend to have rectangular edge geometries, and also rectangular 
shaped heater regions. The rectangular case would be thus interesting to study, but it 
is unfortunately more complex, resulting in fully two-dimensional integral equation. Our 
results should, nevertheless, describe fairly accurately even the rectangular geometry, if 
temperatures are not probed too near the heater or too near the membrane edge. Another 
future extension of this work could also be the inclusion of specular scattering probability 
component, as it is known that for many smooth samples a large portion of the scattering 
on the surfaces is specular—. By increasing the specular scattering for thin membranes, one 
then crosses over to the regime where phonon modes are modified due to interference, and 
become the so called Lamb modes with non-linear dispersion relations^. It is known that 
the Lamb modes affect heat conduction strongly at least in the ballistic limits, therefore 
the cross-over from Casimir to ballistic conduction should also be studied further. 
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